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Abstract 

In this paper, the problem of Magnetic Resonance (MR) image reconstruction from partial Fourier 
samples has been considered. To this aim, we leverage the evidence that MR images are sparser than their 
zero-filled reconstructed ones from incomplete Fourier samples. This information can be used to define 
an optimization problem which searches for the sparsest possible image conforming with the available 
Fourier samples. We solve the resulting problem using the well-known Alternating Direction Method of 
Multipliers (ADMM). Unlike most existing methods that work with small over-lapping image patches, 
the proposed algorithm considers the whole image without dividing it into small blocks. Experimental 
results prominently confirm its promising performance and advantages over the existing methods. 

Index Terms 

Magnetic Resonance Imaging (MRI), image reconstruction, partial Fourier samples, ADMM. 

I. INTRODUCTION 

The pervasive impact of the Magnetic Resonance Imaging (MRI) as a competent approach to investigate 
the anatomy and physiology of the body can not be ignored in the present context of medical imaging. 
In MRI, the data are acquired as partial Fourier samples in k-space. The goal is to accurately reconstruct 
the medical images from these incomplete samples jointly assuring for the minimal scan time. 

Many techniques have already been proposed to reconstruct images using as fewer samples as possible, 
and hence reducing the scan time JT|, Q. One such state of the art approach is to use the theory of 
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Fig. 1. Histograms of original (top) and zero-filled reconstructed (bottom). 


compressed sensing (CS) which demonstrates highly promising results 0, 0- For this, the overall 
reconstruction is formulated as a sparse coding problem, in which the original image is assumed to 
have a sparse representation in an appropriate domain like wavelets, Discrete Cosine Transform (DCT), 
contourlets, etc., while conforming with the available compressed measurements (i.e., partial Fourier 
samples). The basic problem formulation is as follows 

min H'Fxllo subject to F u x = y, (1) 

X 

where x is the vectorized version of the estimate of the original image, ||.||o denotes the to pseudo 
norm which counts the number of nonzero elements, \F is some sparsifying basis like DCT, F. u is 
an undersampled Fourier matrix, and y is the corresponding available measurements (Fourier samples). 
Numerous algorithms exist which aim to solve this problem or its noise-aware variant and get approximate 
solutions, see e.g., Iterative Method with Adaptive Thresholding (IMAT) 0, 0, Q, |8j, Smoothed LO 
norm (SLO) 0, Orthogonal Matching Pursuit (OMP) iflOll . 

In DU, the authors proposed to use adaptive sparsifying transforms instead of fixed ones like DCT 
and wavelets. The idea of dictionary learning in sparse signal processing |[l2ll . lfl3l . lfT4ll is adopted to 
simultaneously learn the sparsifying transform and reconstruct the image. The experimental results show 
significant improvements over previous algorithms. 

In this paper, we propose a simple yet fast and efficient algorithm for the reconstruction of MRI images. 
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The idea is based on the evidence that MR images are sparser than their zero-tilled reconstructed ones, 
i.e., the ones obtained as the inverse Fourier transform of the available Fourier samples by replacing the 
unknown samples with zeroes. This is illustrated in Fig. 1. In this figure, the histograms of one typical MR 
image and that of its zero-tilled reconstruction arc shown. As can be seen, the original image is sparser 
than its zero-tilled reconstruction. This motivates us to formulate the image reconstruction problem as 
looking for the sparsest possible image conforming with the available Fourier samples. We solve it using 
the well-known Alternating Direction Method of Multipliers (ADMM) method. We work with the whole 
image instead of small overlapping blocks as in ifTTll . Our experimental simulations indicate improvements 
in both running time and quality of the reconstructed images, as compared to lUTl . 

The rest of the paper is organized as follows. In Section[TI]we formulate the problem. This is followed by 


a review on ADMM and applying it on our own problem. Simulation results are presented in Section III 


Finally, concluding remarks are presented in Section IV 


II. PROBLEM FORMULATION 

Consider an image Xo £ M. NxN . There are only a subset of its Fourier samples indexed in 0 where 
|Q| = M with M usually being much smaller than the whole samples, i.e., N x N = N 2 . The objective 
is to reconstruct the original image Xo from its incomplete Fourier samples in 0. Hence, we exploit the 
evidence which was discussed in previous section and illustrated in Fig. 1. To measure the sparsity, we 
use the t\ norm ||X||i, which is defined as the sum of the absolute values of the entries in X. 

Consider the Discrete Fourier Transform (DFT) matrix D £ C NxN . From the separability of DFT, 
the 2D DFT of an image X is defined as Y = DXD U , where “H” denotes conjugate transpose^ The 
inverse DFT is simply X = D f/ YD. Using this, we propose the following problem to reconstruct the 
original image Xo (Xo = IFFT(Y)) 

nfin ||D h YD||i subject to Y(fi) = Y 0 (D). (2) 

Note that the noise-aware variant of the above problem is as follows, which we do not discuss in this 
paper, 

mm ||D h YD||i + A||Y(fi) - Y 0 (D)|||. (3) 

To solve ([2]), we use the well-known Alternating Direction Method of Multipliers (ADMM). In priority, 
we first review the main points of ADMM. 

'Note that there is a normalization factor of M in D H which we have omitted for simplicity. 
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A. ADMM 

Consider the following equality-constrained problem 

min/i(U) + / 2 (V) subject to AU + BV = C, (4) 

U, V 

where U and V are optimization variables (not necessarily of the same size), f\ and / 2 are convex 
functions, and A, B, and C are the problem parameters. A common approach to solve such problems is 
to use the notion of augmented Lagrangian function, which is defined as follows 

MU, V, A) = f\ (U) + / 2 (V) - trace(A t (AU + BV - C)) + ||| AU + BV - C\\ 2 F , (5) 

where, A is the Lagrangian multipliers matrix associated to the equality constraint and ^ is a positive 
parameter. The idea of an ADMM solver is to iteratively minimize the above augmented Lagrangian 
function over U and V, along with update of A. 

B. Algorithm 

To solve ([2]) using ADMM, we reformulate it as follows 

min ||Z||i subject to Y(f2) = Yo(O), Z = D^YD, (6) 

Z, Y 

where Z is an axillary variable of the same size as Y. The Lagrangian function for this problem is 

MZ,Y, Ai, A 2 ) = ||Z||i - trace(Af (Y(ft) - Y 0 (D)) - trace(A^(Z - D^YD)) 

+ ^HY(fi) - Y 0 (fi))|||. + f IIZ - D W YD|||, (7) 

The Z-update problem, after some straightforward calculations becomes 

min ||Z||i + ^||z - B h YB - —A 2 |||, (8) 

z 2 /j 2 

This problem is separable over entries of Z, which ends up with N 2 scalar sub-problems of the form 

min |x| -f- — F (x — a) 2 (9) 

x 2 

where x and a denote typical entries of Z and D /, YD + ^ A 2 , respectively. Problem (J9|) is well-known 
in the compressed sensing community, which has the following simple closed-form solution 

x*=S^(a), (10) 

where S\ (a) is the soft-thresholding function defined as 

S\(a) = sign(a) • max(|a| — A, 0). (11) 
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Final update formula for Z is 


Z = 5 ± (D fl YD + iA 2 ), 
« fl2 


( 12 ) 


in which, the soft-thresholding function acts component wise. 

To update Y, after re-arranging the terms in Q, we have the following problem 


min ^||Y(n) - Y 0 (ft) - — A x ||^ + ^||Z - D"YD - — \ 2 f F . 
Y 2 fi\ 2 /i 2 


1 


F2, 


1 


(13) 


We use the orthonormal property of the DFT matrix, i.e., D ;/ D = I, to rewrite problem ( fT3| ) into the 
following equivalent form 


min ^||Y(n) - Y 0 (fi) - —Aif F + ^||Y - D(Z - —A 2 )D / NJ? . 

Y 2 H\ 2 H2 

Setting the gradient of the above objective function with respect to Y equal to 0 results in 

IH P. * (Y(fi) - Y 0 (ft) - —Aj) + a» 2 (X - D(Z - —A 2 )D^) = 0, 

Ft F2 


1 


F2, 


1 


\H 112 


(14) 


(15) 


where P is a binary matrix having 1 in the entries in fl and 0 elsewhere and .* denotes component wise 
multiplication. 


Let us define 


1 


A = D(Z-A 2 )D", 

F2 


and 


B = 


Fi(Yo + ^At) + /x 2 A 
fit + F2 


From ( [15] ), we have 

B (i,j) <G 
A (i,j) 

The formula for updating the Lagrangian multipliers matrices are as follows 


Y (i,j) = 


(16) 


(17) 


(18) 


Ar = A, — fii(Y(Ci) — Y 0 (f2)), 

(19) 

A 2 = A 2 -/i 2 (Z-D w YD), 

(20) 
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Algorithm 1 Proposed algorithm. 
Require: Yo(fi), p\ and p 2 - 

Initialization: Ai = A 2 = 0 


for iter = 1 , 2 , ... do 

Z = Si(IFFT(Y) + ^A 2 ) 

Update Y using ( fl8| ) along with ( |2Tj ) and (T7| ). 

A 1 = A 1 -/r 1 (Y(U)-Y 0 (U)) 

A 2 = A 2 — /i 2 (Z — IFFT(Y)) 


end for 

Output: X = IFFT(Y) 


C. Fast implementation 


Note that the matrix multiplications DZD tf and D ;/ YD can be efficiently implemented using the 
Fast Fourier Transform (FFT) of Z and the inverse FFT (IFFT) of Y, respectively. For example, instead 
of ( |T6l >, we write 

A = FFT(Z-—A 2 ), (21) 

F 2 


and instead of (20) we write 


A 2 = A 2 -/i 2 (Z-IFFT(Y)), 


( 22 ) 


The main steps of the final proposed algorithm are summarized in Alg. [T] 


III. EXPERIMENTS 

In this section, we evaluate the performance of the proposed algorithm. Since a comparison to all other 
methods is not possible, so we chose the method in ifTTl as an example. MATLAB R2014a on a system 
with an Intel Core i5, 2.27 GFIz CPU, 4 GB memory and 68-bit windows 8 was employed to run the 
codes. The MATLAB codes, reference image and masks used in lfTH arc available at the second author’s 
web page0 

For comparison, both the algorithms (i.e., the DLMRI and our proposed one) were subjected to two 
different datasets, i.e., the test MR image used in IfTTl and another MR image, which we name as MRI- 
1 and MRI-2, respectively. These data are depicted in Fig. 2. Various sampling masks i.e., Random, 

^ http://www.ifp.illinois.edu/~yoram/DLMRI-Lab/Documentation.html 
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Fig. 2. Sample results of the algorithms. There are three separate sections, i.e., (a)-(c). In each section, from left to right 
and top to bottom the images correspond to: the used mask, the zero-filled image, output of the proposed algorithm, output of 
the DLMRI algorithm, reconstruction error magnitude for our algorithm, and reconstruction error magnitude for DLMRI. (a) 
Random (5%), the vivo image, (b) Radial (16%), the vivo image, and (c) Random (25%), the vivo image. 


Cartesian and pseudo radial available in ifTTTl were applied on the reference data in k-space domain to 
reach different levels of undersamplings. It is noticeable that in pseudo radial mask, a 512 x 512 Cartesian 
grid is used to map the samples on the radial lines to their nearest points on the grid. 

We considered the time it takes to attain convergence in both the algorithms. To evaluate reconstruction 
qualities, as in ifTITl we used Peak Signal to Noise Ratio (PSNR) in dB which is computed as the ratio 
of the peak intensity value of the reference image to the root mean square (RMS) reconstruction error 
relative to the reference image. 

The parameters of DLMRI were set according to ifTTIl . For the proposed algorithm, we fixed the 
parameter /jj = 10, as it hardly affects the results significantly. The parameter //2 is also scaled in the 
range of 10 to 30, so as to opt the best value for each mask. Subsequently, one of the advantages of our 
algorithm over DLMRI is its fewer number of parameters. 

Table I shows the results on both test images. For comparison, the run time, start up PSNR, i.e., 
PSNR(init.), and its value after reconstruction, i.e., PSNR(end) have been reported. It can be seen that 
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TABLE I 


PSNR IN dB and running times. In each cell, top: our algorithm and bottom: DLMRI fin . 


Data 

Mask 

Vivo MR image [6] 


32-channel k-space data 

type 

percent 

(%i 

Time 

(Minute) 

PSNR 

(end) 

PSNR 

(int) 


Time 

(Minute) 

PSNR 

(end) 

PSNR 

(int) 

Random 

5 

0.47 

35.60 

22.03 


0.05 

38.65 

24.96 

75 

37.59 

110 

37.67 

12.5 

0.17 

42.89 

25.75 


0.05 

44.47 

29.13 

77 

40.82 

88 

44.12 

16.67 

0.17 

44.44 

26.03 


0.05 

44.96 

29.87 

78 

41.56 

71 

44.80 

25 

0.17 

47.14 

29.98 


0.05 

49.10 

33.16 

79 

43.46 

72 

47.56 

40 

0.17 

51.03 

34.38 


0.05 

53.70 

39.81 

80 

46.19 

80 

53.01 

Cartesian 

19.14 

0.17 

40.44 

30.53 


0.05 

42.23 

33.74 

72 

41.31 

71 

42.30 

25 

0.17 

42.84 

29.61 


0 05 

43.72 

33.02 

74 

42.53 

71 

44.30 

39.26 

0.23 

46.53 

42.06 


0 05 

45.94 

43.07 

75 

45.36 

72 

47.29 

Radial 

13.91 

0.17 

36.20 

30 


0.05 

41.65 

33.56 

75 

37.93 

71 

41.55 

16.41 

0.17 

38.28 

29.19 


0.05 

43.16 

35.28 

78 

39.16 

71 

42.74 


our algorithm is significantly faster than DLMRI. Moreover, it outperforms DLMRI by 2-5 dB, mostly 
in case of random sampling. 

Figure 2 shows the sampler mask, zero-filled image, reconstructed image using DLMRI and our 
algorithm, and the reconstruction error magnitude images for some sample situations. It can be seen 
that DLMRI smooths out the images while our proposed algorithm preserves the details. The PSNR 
versus iterations are also demonstrated in Fig. 4. This figure shows that the convergence rate of our 
algorithm is higher than that of DLMRI. 

IV. CONCLUSION AND FUTURE WORK 

In this paper, we proposed a new reconstruction algorithm for MR images from their partial Fourier 
samples, which is governed by the sparsity of the original image. We formulate it as a minimization 
problem and solved it using ADMM. The resulting algorithm has a fast implementation using FFT. Also, 
its another important advantage over state of the art methods like DLMRI ifTTil is its fewer free parameters. 
Experimental results further confirm our arguments. 

For future work, extension of the algorithm to noisy case, improving the solutions using other image 
priors and regularizations like total variation, and applying adaptive sparsifying transforms are of interest. 
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Fig. 3. Test data, (a) the 32-channel k-space, (b) the vivo. 



Fig. 4. PSNR (dB) versus iteration number. 
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